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ABSTRACT 

This work presents an implementation of the resistive MHD equations for a generic algebraic 
Ohm's law which includes the effects of finite resistivity within full General Relativity. The 
implementation naturally accounts for magnetic-field-induced anisotropics and, by adopting 
a phenomenological current, is able to accurately describe electromagnetic fields in the star 
and in its magnetosphere. We illustrate the application of this approach in interesting systems 
with astrophysical implications; the aligned rotator solution and the collapse of a magnetized 
rotating neutron star to a black hole. 

Key words: MHD - plasmas - gravitation - methods: numerical 



1 INTRODUCTION 

Magnetic fields play an important role in the dynamics of many rel- 
ativistic astrophysical systems such as pulsars, magnetars, gamma- 
ray burst (GRBs) and active galactic nuclei (AGNs). In many of 
these scenarios, the Ohmic diffusion timescales of the magnetized 
plasma is much longer than the characteristic dynamical timescale 
of the system, so one can formally take the limit of infinite elec- 
trical conductivity. This is regarded as the ideal MHD limit, and it 
is in general a good approximation to describe astrophysical plas- 
mas. Furthermore, such a limit is described by a relatively manage- 
able, but certainly involved hyperbolic system of equations without 
stiff terms which facilitates its computational implementation. The 
ideal MHD limit has been extensively used in the last years to study 
many of the previous systems (i.e., which basically consist of mag- 
netized neutron stars and black hole accretion disks) in the fully 
non-linear regime. 

In spite of its success and convenience, the ideal MHD approx- 
imation also has some limitations. At a purely theoretical level, the 
assumption of vanishing electrical resistivity prevents some impor- 
tant physical phenomena such as dissipation and reconnection of 
the magnetic field lines. Reconnection efficiently converts magnetic 
energy into heat and kinetical energy in very short timescales. This 
process is believed to be the mechanism originating many energetic 
emissions, such as in soft gamma-ray repeaters (which could be ex- 
plained by giant magnetar flares), the Y-point of pulsar magnesto- 
sphere or even the short Gamma-Ray Bursts fUzde nsky|201 1^ . In 
order to describe such processes, schemes going beyond the ideal 
MHD limit are required. 

At the numerical level, all numerical schemes inherit some nu- 
merical resistivity which depends strongly on the resolution, mak- 
ing difficult to disentangle physical phenomena from numerical 
artifacts especially in highly demanding computational scenarios. 
The presence of magnetic fields demands relatively high resolution 
to accurately capture all the physical processes involved, many of 



them occurring at very small scales. This high resolution is par- 
ticularly important in the case of instabilities which amplify the 
magnetic field, such as the Kelvin-Helmholtz instability occurring 
during the merger of binary neutron stars ( [Price & Rosswog|2006[ 
|Obergaulinger et al.|20I0[ l, and the Magneto-Rotational Instability 
(MRI) occurring in accretion disks ([B albus & Hawley 1991 ; Haw^ 
[ley & Balbus 1991[ [Hawley et al . 1995 Balbus & Hawley 1998]^ 
Accurate modeling of the rarefied magnetospheres of compact ob- 
jects similarly requires high resolution. The electromagnetic fields 
in this region may be easier to model by adopting a different limit 
of the MHD equations known as \he force-free limit ^Goldreich &| 
[Julian 19691. In this approximation the fluid inertia is neglected, 
implying that the fluid does not influence directly the dynamics of 
the electromagnetic fields. 

One possibility to overcome these limitations is to consider 
instead the resistive MHD framework and solve the full Maxwell 
and hydrodynamic equations. The coupling between these two is 
provided by the current -by a suitable Ohm's law-. With a conve- 
nient choice of current, including both induction and Ohmic terms, 
it is possible to recover both the ideal MHD limit, as well as the 
finite-resistivity scheme required to describe physical dissipation 
and reconnections. The effect of small-scales-dynamics can also 
be modeled with moderate resolutions by using a suitable current. 
Finally, magnetically dominated magnetospheres can be described 
by a phenomenological current that decouples the fluid from the 
force-free EM fields. 

The numerical evolution of this resistive MHD code is not free 
of difficulties. The resistive MHD equations can be regarded as an 
hyperbolic system with relaxation terms that become stiff for some 
limits of the current. Consequently, numerical evolution of this sys- 
tem represents a numerical challenge, and several works have re- 
cently explored different possibilities to implement it (Komissarov 
20071 Palenzuel a et al.|2009| [Dumbser & Zano tti 2009, Zenitam 
et al.||2010t ^Takamoto & Inoue||20II| [Bucciantini & Del Zanna 
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|2012[[Dlonysopoulou et al.|2012^ . In this work we take a step fur- 
ther in the development of one of these approaches, based on the 
Implicit-Explicit (IMEX) Runge-Kutta. Our aim is to model both 
the interior and the exterior of a star with a phenomenological cur- 
rent based on physical arguments. This will be particularly inter- 
esting to study the electromagnetic emissions of astrophysical rel- 
ativistic systems involving magnetized neutron stars. 

The capabilities of our approach are tested by consider- 
ing the force-free aligned rotator solution, a well studied prob- 
lem in the context of pulsar magnetospheres (Contopoulos &' 
Spitkovskyl|2006 [ |Spitk ovsky|2006,, McKinney_2006. Bucciantini| 



et al.|2006 



Tchekhovs 



Kala pothara kos & Contopoulos|2009[ |Li et al.|2012 
coy & Spitkovsky||2012^ . These works were restricted 



to flat spacetime and excluded the interior of the star from the com- 
putational domain, thus side-stepped the stiffness problem men- 
tioned above. Recently a hybrid approach, matching both the ideal 
and force-free system of equations, revisited this problem within a 
framework capable of studying both star and surrounding magne- 
tosphere within General Relativity ^ Lehner et al.|20lT) . However 
such a scheme still relies on two different approximations applied 
in two regions. The approach we present here finally allows for 
treating the system from a global point of view with a single, gen- 
eral relativistic, framework. 

We also revisit another important astrophysical scenario with 
a much less understood dynamics; the collapse of a magnetized 
neutron star to a black hole. This system represents even a more 
challenging problem because of the strong gravity fields, and it 
has been studied numerically by considering different approxima- 
tions. An early study matched an analytical solution for the star 
to an electrovacuum magnetosphere ^Baumgarte & Shapiro|2003] l. 
More recently, further realism was achieved by adopting the hybrid 
scheme that matched the numerical solution of the star to force- 
free magnetosphere l |Lehner et al.||2011[ l. A step towards study- 
ing this system within a common, resistive, framework was pre- 
sented in jPionysopoulou et al.|2012} , although the star's exterior 
was treated as an electrovacuum magnetosphere. Our approach pre- 
sented here is able to consistently study the star and its force-free 
magnetosphere within the general relativistic resistive MHD equa- 
tions. 

The paper is organized as follows. Section |2] summarizes the 
fully relativistic resistive MHD system, which is slightly differ- 
ent from the one adopted in ( Dionysopoulou et al.|2012i. In Sec- 
tion[3]it is discussed a generic family of algebraic Ohm's law, and 
how to construct a phenomenological current to recover both the 
ideal MHD and the force-free limits. Section|4] summarizes briefly 
the IMEX Runge-Kutta methods and different techniques to solve 
generically the implicit step for any algebraic form of the relax- 
ation terms. The application of these methods to the resistive MHD 
system is performed is section|5] Section|6]presents our numerical 
results for the aligned rotator and the collapse of a neutron star to a 
black hole. We conclude with some remarks. 

Throughout this work we adopt geometric units such that 
G = c = 1, and the convention where greek indices /i, v, a, ... de- 
note spacetime components (ie, from to 3), while roman indices 
i,j,k,... denote spatial ones. Bold letters will represent vectors. 



2 THE EVOLUTION EQUATIONS 

This section summarizes the general relativistic resistive mag- 
netohydrodynamic equations, that will allow us to model self- 
gravitating magnetized fluids. The evolution of the spacetime ge- 



ometry is governed by Einstein equations. The electromagnetic 
fields and the fluid obey, respectively, the Maxwell and the General 
Relativistic Hydrodynamic equations. The closure of the system is 
given by two constitutive equations; the first one is the equation 
of state, which relates the pressure to the other fluid variables. The 
second one is Ohm's law, -defining the coupling between the fluid 
and the electromagnetic fields- which will be described in the next 
section. 



2.1 Einstein Equations 

The geometry of the spacetime can be obtained by solving the four- 
dimensional Einstein equations. These equations can be recast as a 
standard initial value problem by splitting explicitly the time and 
the space coordinates through a 3-1-1 decomposition, such that the 
line element can be expressed as 

,2 



ds' 



-a-dt^ + Yij (dx' +p'dt)( dx-i + j8 ^ dt 



(1) 



where g^y is the spacetime metric, yij = gjj is the intrinsic metric of 
the spacelike hypersurfaces, and the lapse function a and the shift 
vector j3' relates how the coordinates change between neighboring 
hypersurfaces. The normal to the hypersurfaces is given explicitly 
by 



„^ = 1(1,-/5') 



iif, = (-a,0) 



(2) 



Indices on spacetime quantities are raised and lowered with the 4- 
metric and its inverse, while the 3-metric and its inverse are used to 
raise and lower indices on spatial quantities. 

The rate of change of the intrinsic curvature from one hyper- 
surface to another is given by the extrinsic curvature 



1 

2a 



(3) 



where is the Lie derivative along the vector j3'. 

At any given time, the spacetime geometry is then fully 
defined by the 3 + 1 variables {a.p' ,Yij,Kij}. We adopt the 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of 
Einstein's equations to evolve a suitable combination of these 
fields, in a form very close to the presented in ( [Campanelli et al.] 
[20061 ). 



2.2 Maxwell equations 

The electromagnetic fields follow Maxwell equations, that in their 
extended version can be written as jPalenzuela et al.|2010c^ 



Kn 



(4) 
(5) 



where {_F^^,*F^^} are the Maxwell and the Faraday tensors, 
is the electric current and {(j), i//} are scalars introduced to control 
dynamically the constraints by exponentially damping them in a 
characteristic time 1 / K ( [Dedner et al.|2002^ . When both the electric 
and magnetic susceptibility of the medium vanish, like in vacuum 
or in a highly ionized plasma, the Faraday tensor is simply the dual 
of the Maxwell one. 



(6) 
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where e^^'^P is the Levi-Civita pseudotensor of the spacetime, re- 
lated to the 4-indices Levi-Civita symbol rj^^"l^ by 

1 



livafi 



(7) 



In this case, both tensors can be decomposed in terms of the electric 
and magnetic fields, 

pnv ^ n^'E'' -n''E>'+E>'''"^ Banji (8) 
*F^'' = ni'B' -n^B^ -e^"''^ Ean^ (9) 

such that E^ and B^ are the electric and magnetic fields measured 
by a normal observer . Both fields are purely spatial, that is. 



E^nn=B^i 



:0. 



The covariant Maxwell equations l|4| |3jl can be written, by 
performing the 3-1-1 decomposition, in term of the electromagnetic 
fields and the divergence-cleaning scalars ( |Palenzuela et al.|2010c^ 

as. 



(d,~^^)E' - 
{d,-J^p)B' + 



£'j''Vj{aBk) + aY-'Vj\l/ 
atrKE' - af 
a'VjE' = aq — aKXff 
e'j''Vj{aEk) + ay'jVj^ 
atrKB' 

aV,S' = —axd . 



(10) 

(11) 
(12) 

(13) 



where e'-'* = e'-'^"«ce = /\/y is the three-dimensional Levi- 
Civita pseudotensor. Since F^*^ is antisymmetric, the four- 
divergence of equation ^ leads to an additional equation for the 
current conservation of Maxwell solutions. 



:0. 



(14) 



The electric current can be decomposed into components along 
and perpendicular to the vector n'^ , 



1+r 



(15) 



where q and are the charge density and the current as observed 
by a normal observer . Again, is purely spatial, so J^ny — 
0. The current conservation ^14\ can be expressed, with the 3-1-1 
decomposition, as 



{d,-^p)q + Vi{aJ') = atrKq 



(16) 



Only a prescription for the spatial components J', which will deter- 
mine the coupling between the EM fields and the fluid, is required 
to complete the system of Maxwell equations. This relation, com- 
monly known as Ohm's law, will be discussed in detail in section|3] 

2.3 Hydrodynamic equations 

A perfect fluid minimally coupled to an electromagnetic field is 
described by the total stress-energy tensor 

T^v = [p(l + e) + p]"fi"v+PgMv 



(17) 



where a factor 1 / v 4;r has been absorbed in the definition of the 
electromagnetic fields. Here p is the rest mass density, £ the in- 
ternal energy and p is the pressure, given by a closure relation 
p = p{p,£) commonly known as the equation of state (EoS). These 
fluid quantities are measured in the rest frame of the fluid ele- 
ment. However, to describe the system is usually more convenient 



to adopt an Eulerian perspective where coordinates are not tied to 
the flow of the fluid. The four-velocity describes how the fluid 
moves with respect to the Eulerian observers, and can be decom- 
posed into space and time components. 



(18) 



where corresponds to the familiar three-dimensional velocities 
as measured by Eulerian observers (i.e., v'^n^ = 0). The time com- 
ponent is defined by the normalization relation m'^m^ = —1, such 
that 



W = -n^u'' = {l-Viv')-^/^ , (19) 



where we can now recognize W as the Lorentz factor. 

In summary, the magnetized fluid is described by the phys- 
ical fields (i.e., the fluid variables and the electromagnetic fields) 
plus the divergence cleaning scalars, which form the set of primi- 
tive variables {p,£, p,v' ,E' ,B' ,q,<j),\if) The matter evolution must 
comply with the conservation of the total stress-energy tensor 



Vyr''^ =0, 



(20) 



which can be expressed as a system of conservation laws for the 
energy density U and the momentum density Si, defined from the 
projections of the stress-energy tensor 



U = n^nvT^^ , Si = -n^T^i , Sij = Tij 



(21) 



In addition to the conservation of energy and momentum, the fluid 
usually also conserves the total number of particles. 



(22) 



where pu^' is the baryon number density. This equation is just the 
relativistic generalization of the conservation of mass. 

As mentioned above, it is necessary to specify the EOS to de- 
fine the pressure and complete the system of hydrodynamic equa- 
tions. Along this paper we will consider either the polytropic EoS 
p = Kp^, which is a good approximation to describe cold stars, and 
the ideal gas EoS p = (T — l)pe, which allows for shock heating 
in the fluid. 



2.4 Resistive MHD system 

The evolution of the electromagnetic fields follows the Maxwell 
equations and the conservation of charge, while the fluid fields are 
governed by the conservation of the total energy, momentum and 
baryonic number. In order to capture accurately the weak solutions 
of these non-linear equations in presence of shocks it is important 
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to express them as a set of local conservation laws, namely 
d,{^B') + 4[^(-i3'^B' + + (23) 

<9,(Vr£'') + 4[Vr(-i3'£'-a(e'*^S,--yV))] (24) 



^/(v/r-^-) + 5t[v/y(-/3*0 + aB'^)] (25) 

d,{Vf¥) + dtiVYi-p'^W+aE'')] (26) 

= ^[—a\i/trK + E''{di^a) + aq—aK\if] 

dtiVfq) + ^t[^(-/3*? + «/)]=0 (27) 

5,(Vr£») + dklVf{-p'' + CCv'')D]=0 (28) 



(29) 



= ^[a5'-'/i:,7-5-'(9yce] (30) 
<9,(^5;) + dt[^{-p^Si + aS''i)] (31) 



where we have defined 
D = pW , 

t = hW^ - p + ^{E^- + b'^) - pW . 

Si = hWh'i + EijuE^B'' , 
Sij = hW^ViVj + Yijp 

^E,Ej-BiBj + ^Yij{E^+B^) , 

and the enthalpy h = p{l + e) + p. This form of the relativistic 
resistive MHD equations is basically the same presented already 
in ( [Dionysopoulou et al.||2012|. Another similar formula tion has 



(32) 

(33) 

(34) 
(35) 



also been derived recently ( Bucciantini & Del Zanna|2012^ . Notice 
also that the energy conservation has been expressed in terms of the 
quantity T = U — D to recover the Newtonian limit of the energy 
density. 



3 COUPLING BETWEEN THE EM FIELDS AND THE 
FLUID 

Maxwell and hydrodynamic equations are coupled by means of the 
current , whose explicit form generically depends on the electro- 
magnetic fields and the local fluid properties measured in the co- 
moving frame. Consequently, it is convenient to introduce the elec- 
tric and magnetic fields measured by an observer comoving with 
the fluid, namely = F^^iiv and b^' = *F^^uv. Notice that, since 
£^Ufi = b^Ufi = 0, there are only three independent components. 
The Maxwell and Faraday tensors can therefore be expressed as 

/tMv = u^'e''-u\^'+e^'''"f baup (36) 

*F^''' = u^'b''-u''b^-E^''"P eaup (37) 

and the electric current can be decomposed into components along 
and transverse to , 



where j^u^ = and q is the charge density measured by the co- 
moving observer. The relation with the Eulerian quantities (jTsjl can 
be obtained from 



(39) 



Substituting these results into eq. l |38[ l and using the 3-1-1 decompo- 
sition of the four- velocity, one can write the spatial components of 
the current as 

Ii=Ji = (q + fn^)vi + ji . (40) 

Since the charge density follows directly from the current conserva- 
tion \\6) , the prescription for the three-dimensional electrical cur- 
rent Ji is the only missing piece to completely determine Maxwell 
equations. 



3.1 Generalized covariant Ohm's law 

A standard prescription, known as the Ohm's law, is to consider that 
the current is proportional to the Lorentz force acting on a charged 
particle, implying a linear relation between the current and the elec- 
tric field in the comoving frame. A richer variety of physical phe- 
nomena may be described by including also additional terms pro- 
portional to the comoving magnetic field, leading to a generalized 
covariant Ohm's law of the form. 



(41) 



being a^^ the electrical conductivity of the medium ( [Bekenstein &| 
|Oron|1978) and X a parameter related to the covariant generaliza- 
tion of the mean-field dynamo jBucciantini & Del Zanna|2012^ . 

The electrical conductivity can be calculated either in the 
collision-time approximation ( [Bekenstein & Oron|1978j) or in the 
framework of relativistic charged multifluids i Andersson||20 1 2) , 
leading to the same main results. The tensorial conductivity can be 
written as, 

^e^'^l'uabp) {AD 



^b^'b^ 



1+^2^2 

where the coefficients are given by 
^ = l/R = etr/nig , <J = R/{nee) 



(43) 



Here T,. is the collision or relaxation time, rig is the electron density 
and e and are the electron's charge and mass. In the framework 
described in ( | Anders son|20 1 2\ , R is introduced as a proportionality 
constant in the dissipative force between the two components of the 
fluid. It is easy to check that the first term of the conductivity ^42\ 
leads to the well known isotropic scalar case, while the other two 
represent the anisotropics due to the presence of a magnetic field, 
corresponding to the Hall effect. 

In order to compute the closure relation (j40j it is necessary to 
write the general relativistic Ohm's law in terms of fields measured 
by an Eulerian observer. Let us first consider a simplified Ohm's 
law neglecting both the dynamo effects and the last term in the 
tensorial conductivity ( |42^ , 
a 



-e{evb')b^] 



(44) 



as it has also been used in ( [Zanotti & Dumbser||201l) . It was 
pointed out that this current implies an incomplete Hall effect jAn-| 
|dersson|2012[ l, but it will be enough for our later discussion. Within 
these assumptions, and using that the electric and magnetic fields 
in the fluid frame can be written as 



I^'^u^'q + f 



(38) 



(45) 
(46) 
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it is straightforward to obtain tiie contraction 

j^n^ = j^^[e^n^+e{evb')b^n^'] (47) 

= Y^^[-W{E\)-We{E^Bk){B\)] . 

The prescription for the spatial current ^40\ can be now computed, 
leading to 



or, after some simple manipulations. 



Ji = qvi 



where we have introduced the shortcuts 



W 



W 



[Ei + SijkvJB'' - 
[Bi-Eiji.vjE''- 



(48) 

(49) 
(50) 



It is important to recall that this current accounts not only for 
isotropic resistivity but also for some anisotropic effects induced 
by the magnetic fields. 

In the regime of low magnetization (i.e., p/B^ ^ 1) these 
anisotropic effects are expected to be small, implying ^ <C 1. In 
this limit the third term in the current can be neglected, lead- 
ing to the well-known isotropic Ohm's law. The high conductivity 
of the fluid implies that, in order to get a finite current, the electric 
field measured by the comoving observers must vanish 



:0 



(51) 



This is the ideal-MHD condition, which states that the electric field 
is not an independent variable since it can be obtained via a simple 
algebraic relation from the velocity and the magnetic vector fields. 

The anisotropic effects are expected to be important in mag- 
netically dominated fluids (i.e., p/B^ <C 1). In this limit <^ 2> 1, and 
the second term in the current (j48j can be neglected. In highly con- 
ducting fluids a finite current is recovered only if the electric field 
is perpendicular to the magnetic field. 



e^bfi =E'Bi = . 



(52) 



since the initial assumption of magnetically dominated fluid pre- 
vents the trivial solution b' = 0. In the next subsection it will be 
shown that this relation is one of the constraints of the force-free 
approximation. 



3.2 The force-free limit 



The magnetospheres of magnetized neutron stars ( [Goldreich & Ju-| 
|lian||I969} and black holes immersed in externally sourced mag- 
netic fields ( [Blandford & Znajek|I977) are filled with a low-density 
plasma so rarefied that even moderate magnetic fields stresses can 
easily dominate over the pressure gradients. In this regime, the 
main contribution to the stress-energy tensor comes from the elec- 
tromagnetic part, r^v ~ 7)iv • Allowing by Maxwell equations, the 
total conservation of energy and momentum can be written as 



= Vvr^" 



(53) 



The vanishing of the Lorentz force F^^ Iv leads to an approxima- 
tion known as force-free limit, which is valid only for magnetically 
dominated plasmas with negligible inertia. The spatial components 
of the force-free condition l |53| l, after performing the 3-)-l decom- 
position, are 



qE' 



e'^'^JiBi, = 



-{J% 



B' 



E'Bi = 



(55) 



where we have defined v'^ = e^^'^E jBj^/ B^ as the drift velocity. Sev- 
eral options have been proposed to compute the term j'^Bi^, which 
is crucial to provide a completely explicit relation for the cur- 
rent. For instance, a closed formed for the curr ent can be calcu - 
lated by enforcing the constraint i9,(£'S,) = I Gruzinov||2007 1. 
Another option is to evolve Maxwell equations by considering 
only the drift term of the current \55) , and correct the electric 
field afte r each timestep to satisfy the other force-free condition 



E'Bj = (Komissarov 2004 Spitkovsky 12006'). This approximation 



has been used successfully to study numerically pulsar magneto- 
spheres ( |Spitkovsky|2006^ and jets emerging from black holes with 
an externally sourced magnetic field ( [Palenzuela et al.||2010b|a[ 
|Neilsenetal.|2011t . 

The force-free limit can also be achieved by considering 
an effective anisotropic conductivity with a generic form given 
by ^Komissarov|2004[poesta et al.|2"0T2{|Alic et al.|2012| l 



Bk)B' + x(E--B^)E' 



(56) 



(54) 



where CT|| is the (anisotropic) conductivity along the magnetic field 
lines. The additional term proportional to E"^ — B^ is introduced in 
order to enforce the physical constraint \E\ > The remarkably 
close resemblance between the covariant current ([48} and the force- 
free one \5(i) suggests that both of them could lead to the same so- 
lutions for some limit of the conductivities. However, the force-free 
current l |56[ ) attains a particularly interesting feature; due mainly to 
the assumption of negligible fluid inertia, it does not depend on the 
fluid fields. This means that the EM fields are decoupled to the fluid 
variables, an advantage that could be used to model accurately the 
EM fields in regions where the fluid description is not accurate. 

3.3 A current for the ideal MHD and the force-free limits 

The numerical evolution of the ideal MHD equations typically 
fails in low density regions with high magnetization unless suffi- 
cient resolution is available, a situation that arises commonly in 
the magnetospheres. A standard practice to avoid these failures is 
to maintain a density floor (i.e., the so called atmosphere) in re- 
gions of low density to exploit advanced numerical techniques for 
relativistic hydrodynamics. The density in the atmosphere is much 
smaller than that inside the star, so this approach does not affect the 
star's dynamics. However, in the magnetosphere the fluid inertia 
(and pressure) is typically much smaller than that of the electro- 
magnetic field and one generally encounters numerical difficulties. 
These problems are mitigated by increasing the density in the atmo- 
sphere, effectively decreasing the magnetization in the exterior of 
the star. Although these modifications produce an unphysical mod- 
eling of the plasma in the magnetosphere, one could still solve cor- 
rectly Maxwell equations by using a suitable current that decouples 
the electromagnetic fields from the fluid variables. 

As explained earlier, the covariant current l[48j reduces to the 
ideal MHD limit for high isotropic conductivities (i.e., CT — > oo and 
E, 0), while that the force-free constraint E'Bi = is enforced 
for large anisotropic conductivities (i.e., a,i, ^ oo). This suggests 
that the solutions for the EM fields in both limits can be achieved 
just by changing the anisotropic conductivity, independently on the 
plasma magnetization. Although Ohm's law |48]l is quite general, it 
still couples the EM fields to the velocity. In addition, the parameter 
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is not appropriate to model the fast decay of tiie magnetic field 
with the distance to the source. To overcome these difficulties, and 
in part motivated by the strategy introduced in i jLehner et al.|201 
we introduce the following phenomenological current to include 
both the ideal MHD and the force-free limits, 



f = q[{l-H) 



-Hv',] 



(57) 



+ 



1 + ^2 



where // is a function which vanishes whereas the ideal MHD limit 
is valid, and tends to 1 whereas the force-free limit is more appro- 
priate. The anisotropic ratio ^, which can be reinterpreted from the 
definition E,^b^ = f ^, can be conveniently set to be a constant in 
the region where the force-free limit is valid. The physical condi- 
tion 5^ — £^ > is enforced through a new current term propor- 
tional to an anomalous conductivity x, which only appears when- 
ever < . Overdamping of the electric field is avoided by set- 
ting this anomalous conductivity to the characteristic decay time 
X ~ (o'vT'*-^^) which can be estimated from the time evolution 
of B^-E^. 

Let us consider the particular astrophysical scenario of mag- 
netized neutron stars. The large fluid conductivity, both inside and 
outside the star, is modeled by using a large constant CT ~ 10^. The 
anisotropic ratio, which defines the regions described either with 
the ideal MHD or the force-free limits, is defined as = Ha. This 
choice ensures that the interior of the star (i.e., = 0) is dominated 
by a large isotropic conductivity, reducing the system of equations 
to the ideal MHD limit. The exterior of the star (i.e., /f = 1) is dom- 
inated by the anisotropic terms which enforce the force-free condi- 
tion. The kernel function H is defined such that vanishes inside the 
star and its value becomes unity outside. A smooth transition be- 
tween the inner and the outer region is achieved by using 



(58) 



We typically adopt K ^ 0.001/pa,„, and Po 50 — 400par,„, being 
palm the value for the density of the magnetosphere. 



4 HYPERBOLIC SYSTEMS WITH RELAXATION 
TERMS 

The general system of relativistic resistive MHD equations ( |23^ 
|31|57^ brings about a delicate issue when the conductivity in the 
plasma undergoes very large spatial variations. In regions with high 
conductivity, in fact, the system will evolve on timescales which 
are very different from those in the low-conductivity region. Math- 
ematically, therefore, the problem can be regarded as a hyperbolic 
system with relaxation terms which requires special care to cap- 
ture the dynamics in a stable and accurate manner. The prototype 
of these systems can be written as 



d,V = F(U) + -R{V) 



(59) 



where e > is the relaxation time. In the limit e — > oo the system 
is hyperbolic with spectral radius c/, (i.e., the absolute value of the 
maximum eigenvalue). In the other limit £^0 the system is clearly 
stiff since the time scale of the relaxation (or stiff term) R{V) is 
much smaller than the maximum speed c/, of the hyperbolic part 
F(U). 

In the stiff limit (e — > 0) the stability of an explicit time evolu- 
tion scheme is only achieved with a time step size Af ^ e, a much 



stronger restriction than the CFL condition At ^ Ax/c/, of the hy- 
perbolic systems. The development of stable and efficient numeri- 
cal schemes to overcome this restrictive constraint is challenging, 
since in many applications the relaxation time can vary many orders 
of magnitude. 

Different alternatives to deal with the inherent stiffness of the 
relativistic resistive MHD equations has been proposed in the last 
decade; combination of splitting methods and analytical solutions 
l |Komissarov 2007; Ze nitani et al.|20 10; Takamoto & Inoue 201l]l, 
discontinuous Galerkin methods (Zanotti & Dumbser 201 1 , Dumb^ 
|ser & Zanottil[2b09|l and Implic it-Explicit (IMEX) Runge-Kutta 



methods ( [Palenzuela et al.|2009| |Bucciantini & Del Zanna|2012[ 
|Dionysopoulou et al.|2012| . The following subsections summarize 
the IMEX Runge-Kutta schemes, a family of time integrators which 
are able to deal with the potentially stiffness issues and are rela- 
tively easy to incorporate into an existing relativistic ideal MHD 
code. 



4.1 Implicit-Explicit Runge-Kutta methods 

An efficient way to solve the hyperbolic-relaxation systems is based 
on the IMEX Runge-Kutta methods. Within this scheme, all the 
fields are evolved by using a standard explicit time integration ex- 
cept the potentially stiff terms, which are evolved with an implicit 
time discretization. For the generic system l |59[ ) this scheme takes 
the form ( |Pareschi & Russo|200"5| l 



A/Xa-yF(uW) 



n+l 



:\J" 



A/X«.7-^(U<^'') 

!=1 !=1 ^ 



(60) 



where U^'^ are the auxiliary intermediate values of the Runge- 
Kutta. The coefficients can be represented as v x v matrices A = 
(ciij) and A = such that the resulting scheme is explicit in F 
(i.e., fly = for j ^ /) and implicit in R. An IMEX Runge-Kutta is 
characterized by these two matrices and the coefficient vectors ro,- 
and CO,-. Notice that at each substep the auxiliary intermediate values 
U^'^ involves solving an implicit equation. Since the simplicity and 
efficiency of solving the implicit part at each step is of great im- 
portance, it is natural to consider diagonally implicit Runge-Kutta 
(DIRK) schemes (a/y = for j > i) for the stiff terms. A deeper dis- 
cussion on the IMEX schemes and the detailed form of the schemes 
considered here are presented in appendix [A| 

4.2 Solving generic systems with IMEX schemes 

The vector of evolved fields U can be split in two sets of variables 
(V, W), depending on whether or not they contain any relaxation 
term in their evolution equations. The evolution system can then be 
generic ally written as 



(9,W 



Fi/(V,W)- 



:^v(V,W) 



(61) 
(62) 



where we have considered that the relaxation parameter e can be 
any function not depending directly on the present value of the V- 
fields. The procedure to compute each auxiliary step U''' can be 
split in two stages: 
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(i) compute first the intermediate values {V*,W*} which in- 
volves information from previous steps, 

(-1 

W*=W" + At "ijFwiV^-'^) 
j=i 
1-1 

V*=V" + Ar £ ciijFviv'^j^) 
i=i 

(ii) include the relaxation term at the present time by solving the 
implicit equation 



W* 



e(') 

which clearly involves only the V-fields. 



(64) 



The complexity of inverting this implicit equation depends on 
the particular form of the relaxation terms. From now on we will 
restrict ourselves to the algebraic case Rv{V) = /(U). Next it is 
described two different ways to solve this implicit equation; the first 
one can only be applied when Ry (U) is a linear function, whereas 
the second one allows Rv(V) to have any non-linear dependence. 



4.2.1 Rv depending linearly on \ 

The simplest case, however enough to cover a broad range of inter- 
esting situations, is to consider a linear relaxation term 



i?v(V,W) =A(W)V + 5(W) . 

The implicit equation ^64\ can then be trivially solved 



(65) 



M 



M 



V*+a,-,-|^5(W«; 
Ar 



(66) 



The matrix inversion can be performed analytically and written in a 
compact form for most of the interesting cases, so that the implicit 
step can be solved in a completely explicit way. 



4.2.2 Rv depending non-linearly on V 

In the more general case -with an arbitrary non-linear dependence- 
it is usually not feasible to solve analytically the implicit step, re- 
quiring some approximation to find the solution. A convenient ap- 
proach to solve this problem is to linearize the stiff term around an 
approximate solution {V, W'''}, such that 



/;y(vW,w«) 



+ 



dRv 



(67) 



v.ww 



(VW_v) 



Notice that we are linearizing around the solution W''', which is 
already known at the beginning of the implicit step. 

and substituting the previous 



By defining A ^(^^.^,,,, 
expansion \bl) in \64) , it is obtained 

V('' = V* +a,7 %-ARv{y) +A(V(') 



-V)] 



(68) 



This implicit equation can be written, after some manipulations, in 
the following way 



M ^ [/_„, iLA(V,wf'))]-' 

£[') 



(69) 



The final expression \69) can be solved through a Newton-Raphson 
iterative procedure such that, at each iteration m, uses an initial 



guess V = v(^_j, to find the next approximate solution V| 



5 NUMERICAL EVOLUTION OF THE RESISTIVE 
MAGNETOHYDRODYNAMICS SYSTEM 

We adopt finite difference techniques on a regular Cartesian grid 
to solve the problems of interest. To ensure sufficient resolution is 
achieved in an efficient manner we employ adaptive mesh refine- 
ment (AMR) via the HAD computational infrastructur e!^ tha t P™" 
vides distributed, Berger-Oliger style AMR ( |Liebling]|2002^ with 
full sub-cycling in time, together with an improved treatment of ar- 
tificial boundaries ( [Lehner et al.|2006| l. The refinement regions are 
determined using truncation error estimation provided by a shadow 
hierarchy ( [Pretorius|2002^ which adapts dynamically to ensure the 
estimated error is bounded within a pre-specified tolerance. The 
spatial discretization of the geometry is performed using a fourth 
order accurate scheme, while that High Resolution Shock Captur- 
ing methods based on the HLLE flux formula with PPM recon- 
struction are used to discretize the resistive MHD variables jAn-| 
[derson et al.|2006||2008^ . The time-evolution is performed through 
the method of lines using a third order accurate Implicit-Explicit 
Runge-Kutta integration scheme described in the previous section. 
We adopt a Courant parameter of A = 0.25 so that A?/ = 0.25Ax/ 
on each refinement level /. On each level, one therefore ensures 
that the Courant-Friedrichs-Levy (CFL) condition dictated by the 
principal part of the equations is satisfied. 



5.1 Evolution of the electric field 

The relaxation terms of the resistive MHD system are associated 
to the current, which mainly appears in the time evolution equation 
of the electric field. The evolved fields can then be split into and 
non-stiff W = {£), T,5,-,B', ^,^,q} and potentially stiff V = {£'}. 
The evolution of the non-stiff fields is performed by the explicit 
part of the IMEX Runge-Kutta, and it is very similar to a standard 
implementation of the ideal MHD equations. The evolution of the 
electric field contains in addition the relaxation terms, namely 

^,(/rE) = FE + {,/yRE) ■ (70) 

Fe = -5t[Vr(-j8'£'-a(e'*-'S,-yV))l, 

Re = -aJl . 

where the factor 1/e, corresponding to the fluid conductivity, is ab- 
sorbed in the definition of Re . The current has been split into a po- 
tentially stiff part, /j, and the terms which can be treated explicitly. 



publicly available at http://had.liu.edu 
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Jl. For the phenomenological Ohm's law lj57| these components 
can be written exphcitly as 



f 



<7 



-H)v'+Hii] 



(71) 



l + ?2 



i?2 



{{e''b,)b'+x{e'--b^)e'} 



Notice that, although the evolution of q is driven by the current, 
these terms do not become potentially stiff in this equation since 
they are not proportional to the field itself. However, the delicate 
balance between the different fields in the current, which allows 
to get finite values even for very high conductivities, may be bro- 
ken during the reconstruction of the fields at the interfaces. These 
unacceptable large errors are prevented in the standard implemen- 
tations of the force-free equations by computing the charge density 
from the constraint q = ViE' instead of using the charge conser- 
vation. The resulting set of equations is still hyperbolic, since the 
charge density onl y couples to the EM fie lds throughout the non- 
principal term qv' ( Palenzuela et al.|2011 1. Here we prefer to keep 
the charge density as an evolution field and treat all the fields in 
the same manner. The errors at the interfaces are avoided by per- 
forming directly the reconstruction of the current J' , which is com- 
puted just after solving the stiff terms. This ensures that the fluxes 
of q will remain bounded between the values given by well-defined 
neighboring points. 



5.2 Inversion from conserved to primitive variables 

The numerical evolution of the resistive MHD system ( |23|31[ > in- 
volves the recovery, after each timestep, of the primitive fields 
{p, e,p,v' ,E' ,B' ,\if,(j) ,q} from the conserved or evolved fields 
y/Y{0, T,Si,E',B' , Y.^.q}. Although the conserved fields are just 
algebraic relations of the primitive ones, the opposite is not true; 
due to the enthalpy and the Lorentz factor these quantities are re- 
lated by complicated equations that can only be solved numerically, 
except for particularly simple equations of state. 

The solution at time f = (n + 1 ) A/ is directly obtained, for most 
of the conserved quantities, by evolving their (non-stiff) evolution 
equations. However, the explicit evolution of the potentially stiff 
fields only provides a partial solution. As explained in the previ- 
ous section, a complete solution for the electric field involves tak- 
ing into account the relaxation terms by solving the corresponding 
implicit equation. For a generic Ohm's law, these relaxation terms 
will depend on the velocity and other primitive fields. Nevertheless, 
the recovery of the primitive variables from the conserved ones in- 
volves all the fields, including the electric field. This is a consis- 
tency constraint which implies that the recovery process and the 
implicit step evolution must be solved at the same time. We will 
next describe an iterative procedure to evolve the stiff part and re- 
cover the primitive fields for the phenomenological current \51\ , as 
described in subsection l4.2.2l 

(i) To start the iterative process it is required an approximate so- 
lution -initial guess- for the electric field £, and the fluid unknowns 
of the system, that we have chosen to be the single combination 
X = hW^. The initial guess for this unknown is given simply by the 
previous time step x = x'"' . Possible choices for the electric field 
initial guess are: 

in) 

• the previous time step = Ej 

• the ideal MHD limit £; = —Eiji^v-'B'', which involves per- 
forming first the recovery in the ideal MHD case (see appendix 
B for details). 



• the approximate solution given by the explicit and previous 
implicit step evolutions = E*. 

• the trivial case E, = 0. 

It may be difficult to estimate a priori which initial guess is more 
convenient. For this reason, our scheme starts with the first option 
and, if no solution is found, tries sequentially the other choices. 

(ii) Subtract the electromagnetic contributions from the energy 
and momentum densities, 



T-^-{E'Ek+B'Bk) 

ok 



Si — Si — EijkE-'B 
such that the Lorentz factor can be computed as, 



x^ - S'Si 



1 

IV2 



= 1- 



(72) 
(73) 

(74) 



(iii) Write also the pressure as a function of the conserved vari- 
ables and the unknowns. For the ideal gas EOS p = (T — l)pe this 
relation is just 



r- 1 



P- 



X 

V/2 



(75) 



(iv) Obtain an equation /(x) = 0, written in terms of the un- 
known X and the conserved fields, such that it is satisfied only for 
true solutions of x. By using the previous expression l |75[ l in the 
definition of f, we can write 



(^lx+[^ 



ilo- 



(76) 



where W is computed through eq.i 74 >. The equation /(x) = 
can be solved numerically by using an iterative Newton-Raphson 
solver. The solution in the iteration m + 1 can be computed as 



"^(m+I) 



(77) 



The derivative of the function /(x) can be computed analytically. 



rx2 

r-i)c- (r-i)z)s2 



r ^rx3 

(v) Update the primitive fields by using the relations 



(78) 



X 

r- 1 



r 



(h-p) 



-52 
p -- 



X 



D 
W 



(79) 



(vi) Update the electric field -with the updated values of the 
primitive fields- by solving the implicit equation, corresponding 
to eq. (|64[l. 



E' = eI + an At R'l; , 



(80) 



whi ch can be formally solved with the method described in subsec- 



tion 



M 



4.2.2 



for V(') = E', that is. 



= E'+M[E*'-E'+aiiAtR'i. 



[I-aiiA,A]- 



54 
dEJ 



(81) 
(82) 



For the phenomenological Ohm's law l |57| ), the matrix M to be 
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inverted is 



■S'r 



W{S'-v'vj) 



(83) 



+ ^{b'Bj + x[2E'Ej + 5'j{E^-B')]} 



with a = an At a a / {I + ^^). 

(vii) Iterate until the solution {x,E'} satisfies their constitutive 
equations f{x),f{E') ^ 10^"*, being /(-£') defined by equation 

In occasions the recovery procedure is unable to find a physi- 
cal state for a given set of conserved variables. In such cases, which 
usually occur near a star's surface, failures can be avoided by as- 
suming that the fluid is isentropic in that timestep and therefore 
satisfying a polytropic EoS p = Kp^ . Since the internal energy is 
also a function of the density (i.e., pe = p/(r— 1)) for isentropic 
processes, the conserved quantities are overdetermined and the en- 
ergy equation can be neglected in the recovery procedure, leading 
to a more robust algorithm. 

Notice also that, although our discussion was focused on the 
phenomenological the Ohm's law jSTj, the method described in 
subsection |4.2.2| can be applied to any algebraic form of the current. 
Even more general cases with derivative terms can be considered, 
with the condition that those must be evaluated at earlier times. In 
a similar way, the method for linear relaxation terms described in 
subsection |4.2. 1 1 can be generically used for non-linear algebraic 
currents with the condition that the non-linear terms are evaluated 
at previous time steps, as it was considered in ( |Alic et aL[T0l2) . 
This option does not require an initial guess for the electric field 
and therefore may be more effective in avoiding unphysical states. 



6 NUMERICAL SIMULATIONS 

In this section we report our numerical studies of astrophysical sce- 
narios involving the dynamical evolution of a rotating magnetized 
star and its magnetosphere. The initial data of rigidly rotating neu- 
tron stars is provided by the LORENE package Magitar|^ which 
adopts a polytropic equation of state P = Kp^ with F = 2, rescaled 
to = 100. Because the fluid pressure in a neutron star is many 
orders of magnitude larger than the electromagnetic one, moderate 
magnetic fields will have an insignificant effect on both the geom- 
etry and the fluid structure, and so they can be specified freely. For 
this reason we have chosen an initial poloidal magnetic field inside 
the star that becomes dipolar in the external region. The electric 
fields are set by assuming the ideal MHD condition, with an initial 
zero fluid velocity in the magnetosphere. 

During the evolution, which is performed with the methods 
described in the previous sections, the ideal MHD and the force- 
free limits are enforced inside/outside the star by using the phe- 
nomenological current \51\ . We monitor the electromagnetic lumi- 
nosity, constructed from the Newman-Penrose scalar <I>2 l |Newmaii| 
|& Penrose] 1962| >, 

J rem » 

iem = --r- = l™ / A'^2?-d£l . (84) 

at r->ooj 

that accounts for the energy carried off by outgoing waves to in- 
finity and it is equivalent to the Poynting luminosity at large dis- 
tances. Additionally we monitor the ratio of particular components 
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Figure 1. Aligned rotator. Several quantities displayed in the equatorial 
plane as a function of the cylindrical radius after two rotational periods. 
The kernel function H indicates the value of the density at which the cur- 
rent changes abruptly. The EM quantities do not display any significant 
discontinuity in that region, as can be appreciated for instance in the charge 
density. The magnetic fields in the magnetosphere, up to the light cylinder, 
co-rotates with the frequency of the star Q.fjs- The anomalous resistivity 
appears only in the regions with E > B, close and beyond the light cylinder. 



of the Maxwell tensor £lf = F,i-/F,-ij, which, in the stationary, ax- 
isymmetric case, can be interpreted as the rotation frequency of the 
electromagnetic field jBlandford & Znajek|1977^ . 

6.1 The aligned rotator 

We consider first the evolution of an uniformly rotating sta- 
ble star of mass M = 1.58M0 and equatorial/polar radius R — 
16.1/10.6 km. The star rotates with a period T = 1.3ms, so that 
the light cylinder is located at Ri^c = cj fl/vij = 62 km. The strength 
of the magnetic field at the pole is Bp = 1.8 x lO'^^G. The numer- 
ical domain extends up to L = 300 km and contains four centered 
FMR grids with decreasing sizes (and twice better resolved) such 
that the highest resolution grid has Ar = 0.76 km and extends up to 
76 km (i.e., beyond the light cylinder). 

This initial configuration is evolved until that the solution re- 
laxes to a quasi-stationary state. Different quantities are plotted 
along the equatorial plane in fig.[T|and that both the initial and the 
final magnetic field solutions are displayed in fig.|2] The relaxed fi- 
nal state has the characteristic features observed in previous works. 
The magnetic fields are being dragged by the fluid rotation in the 
interior of the star (i.e., as in the initial state), producing a tension 
that forces the magnetic fields in the magnetosphere to co-rotate 
with the star up to the light cylinder. Beyond this surface, the mag- 
netic field lines open up, creating a current sheet in the equatorial 
plane where the anomalous resistivity in the current (or bringing 
back the neglected fluid inertia) is necessary to preserve the physi- 
cal condition B^ > E^. 

We have computed the Poynting-vector luminosity at two 
surfaces at R^xt = {76,114} km located outside the light cylin- 
der, where the measures converge to a unique well-defined value. 
The EM radiation is mainly dipolar (i.e., around 90% of the en- 
ergy), with a small fraction in higher multipoles. The luminosity 
can be compared with previous results in flat spacetime geometry 
where the spherical star is modeled through inner boundary condi- 
tions l |Contopoulos & Spitkovsky|2006||Spiflcovsky|2006^ 



publicly available at http://www.lorene.obspm.fr 
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Figure 2. Aligned rotator. The fluid density, the magnetic field -poloidal components in lines and toroidal one in colors- and the coefficient x of the anomalous 
conductivity on the x = plane at f = and after two rotational periods of the star The relaxed solution exhibits the known properties of the aligned rotator 
solution, namely an opening of the magnetic field lines roughly at the expected position Rlc ~ 4.0 S,. These plots do not show the entire computational 
domain. 



Our results agrees within a difference of ~ 20%, where we have 
used ^iMs = R^q. It is unclear where this small disagreement may 
come from, since there are several possible explanations; the am- 
biguity in the definition of the radius of oblated stars, an excess of 
dissipation in the current sheet, or purely strong gravitational ef- 
fects, which may become important due to the high compactness 
M/.R = 0.125 of the star. 

We have also monitored both the energy-momentum con- 
straints and the divergence constraints, checking that they re- 
main small and under control during the evolution. In particular, 
|V • 5|/|6| ^ 0.05 in all the domain but the current sheet. By com- 
paring the solutions obtained with three different resolutions, each 
one improving a factor 1.18 the previous space discretization Ax, 
we have observed that the code converges at 1.8-order. The lumi- 
nosity for these three resolutions displayed in fig. [3] shows that, 
in spite of the spasmodic reconnections happening in the current 
sheet, the system converges to a quasi-stationary solution with a 
steady luminosity. 



6.2 Collapse of a magnetized rotating neutron star 

After assessing the validity of our implementation with the aligned 
rotator solution, we can consider a more challenging and dynami- 
cal case; the collapse of an uniformly rotating magnetized neutron 
star to a black hole. The initial data is the same as it was considered 
in ( [Lehner et al.||20lT| l; a star lying on the unstable branch with 
mass M = 1.84M0 and equatorial/polar radius R = 10.6/7.3 km, 
rotating with a period T = 0.78ms so that the light cylinder is lo- 
cated at ^LC = 37 km. The strength of the magnetic field at the 



rescaled to any strength as long as the magnetic pressure is much 
smaller than the fluid one. The numerical domain extends up to 
L = 300 km and contains 6 centered FMR grids with decreasing 
sizes such that the highest resolution grid has Ax = 0.19 km and 
extends up to 21 km, while that the second highest extends up to 
44 km, beyond the initial location of the light cylinder. 

Small perturbations arising from numerical truncation errors 
are enough to trigger the collapse of the unstable star. The horizon 
appears after around Ims, although the most dynamical part only 
stands for the last 0. Ims, ending when all the matter disappears be- 
yond the horizon and the nearby magnetic fields reconnects in the 
equatorial plane and escapes to infinity. The conservation of angu- 
lar momentum implies that the angular velocity of the star increases 
during the collapse, dragging the magnetic field lines in the mag- 
netosphere and bringing the light cylinder closer to the star. The 
magnetic fields also grow due to the magnetic flux conservation. 
Once all the fluid has accreted onto the black hole, the magnetic 
fields looses their anchorage, reconnects and propagates away from 
the source. A significant fraction of the energy stored in the mag- 
netosphere is radiated to infinity in this burst. The density of the 
star, the Poynting vector density |<52p ttid the magnetic fields are 
displayed at some representative stages of the collapse in fig.|4] 

The growth of the angular velocity and the magnetic field im- 
plies that the luminosity of the aligned rotator l |85t during a quasi 
adiabatic collapse will increase as Lq (i^NS |Lyutikov|20I 



being Lq the initial luminosity of the star. However, since the col- 
lapse time is shorter than the star's period, the outer part of the 
magnetosphere is not able to respond to the changes in th e start's 



surface, reducing the power of the luminosity to (R^^s/R) (Lehner 



pole is chosen to be Bp = 1 .8 x lO" G, although the results may be |et al.|20l"T| l. In addition, strong gravitational effects will soften the 
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Figure 3. Aligned rotator. The EM luminosity as a function of tlie rotational period for three different resolutions Ax = {0.76. 0.64, 0.55}km, showing the 
initial transient followed by a fast decay to the quasi-stationaiy solution. The luminosity has been normalized with respect to the asymptotic value, reached 
approximately after 2 rotational periods, of the low resolution simulation. 
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Figure 4. Collapse of a magnetized rotating star. The fluid density, the Poynting flux density and the poloidal magnetic field in the x - 
stages of the collapse corresponding to i = {-0.43,-0.30,-0.18,-0.05.0.08,0.20} ms. 



■ plane representative 



growth of both the angular frequency and the radial magnetic field, 
leading to a much more moderate luminosity growth. 

We have computed the electromagnetic luminosity in a sphere 
located at R^^, = 76 km, beyond the light cylinder. The EM radia- 
tion is mainly dipolar and grows during the collapse, with a strong 
burst due to the reconnection when the fluid is completely swal- 
lowed by the black hole. The luminosity and the angular velocity 
- computed inside and outside the star- are displayed in figure |5] 
The energy in the magnetosphere increases by a factor Cp^aji ~ 2 
during the collapse. The total radiated energy can be expressed as 
a fraction of the peak energy CpeaiiEijipgigX), namely 

2 



Emd ~ 1.4 X 10'*''Cpeak El-ad 



lOi^G 



erg. 



(86) 



where we have used iidipole,0 = 1-4 x lO^^S^g 15 ^''S f*"" ^ ^t^"" 
of radius Rj^s ~ 12 km ( [Lehner et al.|[20lT| l. In our simula- 
tion we have found e^,^^ = 0.6, implying that the system radiates 



^-rad ~ 1-6 X 10^^ ergs during the collapse (for a magnetic field of 
lO'^G). Notice that this value is different from the analytical esti- 
mates and indicates the importance of the fast dynamic and strong 
gravitational effects in this scenario. 



7 SUMMARY 

We have presented a formulation of the general relativistic resis- 
tive MHD equations. We have discussed different generalizations 
of the isotropic Ohm's law, and constructed a phenomenological 
current such that the system reduces either to the ideal MHD limit 
or to the force-free approximation just by changing the ratio of 
isotropic/anisotropic conductivities. We have explained how to deal 
with the potential stiffness of the equations by using the implicit- 
explicit Runge-Kutta methods, showing how to perform the im- 
plicit evolution of the electric field and the recovery of the prim- 
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Figure 5. Collapse of a magnetized rotating star. The EM luminosity and tlie angular frequency of the magnetic field - computed inside and outside the 
star- as a function of time. These quantities have been normalized with respect to some reference values, calculated once the system relaxes to a quasi-steady 
state (i.e., the aligned rotator solution) at early times. The agreement between the interior and the exterior angular velocity shows that it is being propagated 
con'ectly through the surface of the star. 



itive from the conserved fields at the same time for any algebraic 
Ohm's law. We implemented the formulation within the HAD com- 
putational infrastructure and revisited two interesting astrophysical 
problems; the aligned rotator and the collapse of a rotating neutron 
star to a black hole. None of these cases has a known analytical so- 
lution, although the first case has been studied extensively. We find 
a reasonable agreement between our results and previous studies 
of the aligned rotator, recovering the same qualitative features and 
approximately the same electromagnetic luminosity. 

The case of the collapsing star is more challenging and has 
been only studied previously either assuming an electrovacuum 
magnetosphere and/or by matching the exterior to the interior solu- 
tion. Our results are qualitatively similar to those found in ( |Lehner| 
|et al.|20lT| l, although the total radiated energy in our simulations is 
one order of magnitude larger due to an increase in both the peak 
energy in the magnetosphere and the fraction of radiated energy. 
The possible detectability of this burst has been already discussed 
in detail in ^Lehner et al.|20ir| l and therefore will not be repeated 
here. 

In conclusion, the resistive MHD framework allows to con- 
sider a broad range of new phenomena; study reconnections and 
dissipation with more realistic Ohm's law - like the resistive so- 
lutions of pulsar magnetospheres l |Li et al.|2012^ -, model the mag- 
netic growth due to different instabilities by using the mean-field 
dynamo ( Bucciantini & Del Zanna|2012[ l, and compute the magne- 
tosphere interaction of binary systems -like neutron-neutron stars 
and neutron-black hole-, which may be crucial to study the possible 
electromagnetic counterparts to the gravitational waves emitted by 
these systems, among others possibilities. Work on these directions 
is in progress and it will be reported in the near future. 



where the coefficients c and c used for the treatment of non- 
autonomous systems are given by the following relation 

1 ! 

Solutions of conservation equations have some norm that de- 
creases in time. It would be desirable, in order to avoid spurious nu- 
merical oscillations arising near discontinuities of the solution, to 
maintain such property at a discrete level by the numerical method. 
The most commonly used norms are the TV-norm and the infinity 
norm. A scheme is called Strong Stability Preserving (SSP) if main- 
tains a given norm during the evolution ( Spiteri & Ruuth|2002| >. 

In all these schemes the implicit tableau corresponds to an L- 
stable scheme (that is, co^A^'e = 1, being e a vector whose compo- 
nents are all equal to 1), whereas the explicit tableau is SSPfc, where 
k denotes the order of the SSP scheme. We shall use the notation 
SSPA:(i, CT,p), where the triplet {s,a,p) characterizes the number 
of s stages of the implicit scheme, the number (7 of stages of the 
explicit scheme and the order p of the IMEX scheme. 

There are different IMEX RK schemes available in the liter- 
ature. We have considered only third order IMEX schemes, some 
of them found in the literature ( Pareschi & Russo 2005 1 and others 
developed by us. All of them are based on a third order SSP explicit 
scheme that can be implemented efficiently by using only two lev- 
els of fields and one of rhs. It is worth mentioning that these meth- 
ods are still under development and have few drawbacks. Probably 
the most serious one is an accuracy degradation for some range of 
the relaxation time e. 



APPENDIX A: IMEX 

IMEX Runge-Kutta schemes can be represented by a double 
tableau in the usual Butcher notation ( [Butcher] 1 987| [2003l l 

c \ A c \ A 



APPENDIX B: IDEAL MHD LIMIT 

The ideal MHD limit can be obtained by requiring the current to 
be finite even in the limit of infinite isotropic conductivity, leading 
to the condition E' = —e'j^VjB]^. The Ohm's law current becomes 
undetermined (i.e., an infinite conductivity multiplying a vanishing 
electric field in the co-moving frame), but it can still be computed 
from the redundant Maxwell equation for the electric field evolu- 
tion ijTTJ. The evolution of the magnetic field can be simplified by 



© 2012 RAS, MNRAS 000,[T]{T4] 



Modeling magnetized neutron stars using resistive MHD 13 

Table Al. Tableau for the explicit (left) implicit (right) IMEX-SSP3(4,3,3) L-stable scheme 
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a = 0.24169426078821 , /3 = 0.06042356519705 . 77 = 0.12915286960590 



Table A2. Tableau for the explicit (left) implicit (right) IMEX-SSP3(5,3,3) L-stable scheme 
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1/6 
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041 = — (2a^ + 2a-l) , 042= — (-4a^ + l) , a43 = -(-3a+l) . a = l/3 
8a 8a 4 



substituting the ideal MHD condition in \13) , 



(Bl) 



The transformation from conserved to primitive is simphfied 
by eliminating the electric field as an independent variable and may 
allow us to recover the primitive quantities in a more robust way. 
Substituting the ideal MHD condition in the definition of the con- 
served variables 

BK 

(B2) 



hW-+B'--p-D-^-[{B''vkf- 



T 

5; = [hW'-+B^]vi-{B''vk)Bi 
it is easy to check that 



SiB' 
MV2 



(B3) 



(B4) 



Using this relation, the scalar product S'S; can be solved for the 
Lorentz factor, obtaining 

x^S^ + (2x + B^){SiB'f 



1 

iy2 



1- 



x2(x + S2)2 



(B5) 



Assuming an ideal gas EoS, and after some manipulations in 
the definition of T (|B2l, the resulting final equation to solve is 



2x2 



{SiB'f 



(B6) 
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